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ABSTRACT 



Optical scintillation on folded paths has been investi- 
gated for the case of spherical waves and a plane mirror as 
a folding target. A criterion for the statistical indepen- 
dence of the two parts of the folded path was formulated, 
utilizing the path integral technique and the Von Karman 
spectrum. Employing the same approach, an expression 
for the scintillation strength (the log-intensity vari- 
ance) on folded paths has been developed. It has been 
shown that the scintillation strength on folded paths 
can be expressed as a product of the one-way scintilla- 
tion strength and a function of the distance source-detector 
in units of the outer scale. The analytical prediction shows 

that the ratio of the log intensity variances on the folded 

17/6 

paths to the one on the one-way path, goes from 2 ' for 

exact folding (the distance source-detector is zero) , to 

11/6 

2 ' when the detector is far away from the source. The 

theoretical predictions have been corroborated by preliminary 
field experiments. 
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I. 



INTRODUCTION 



Intensity fluctuations of visible or near-visible sources, 
while viewed through the atmosphere (scintillation) , play a 
major role in a large variety of situations such as star 
observations, optical communication channels, Lidars, power 
transmission, etc. It has been understood for a very long 
time that the stars' twinkling is caused by the inhomogeneity 
of the atm.ospheric index of refraction. It is only in the 
last thirty years, and especially since the first lasers were 
built about twenty years ago, that this intuitive understand- 
ing has been put into a rigorous quantitative model. Some 
of the better known contributors to this progress are 
Tatarski [1,2], Fried [3], Lutomirski and Yura [4], Ishimaru 
[5], Lee and Harp [6], Clifford [7], Fante [8,9] and Dashen 
[10] . As a result of their work and many others it is now 
believed that in most of the situations the processes involved 
are fairly well understood and the theoretical predictions 
are reasonably close to experimental findings. 

Among the problems which are less understood and for which 
existing theory either does not give reasonable predictions 
or has not been confirmed experimentally, is the problem of 
scintillation on folded paths. In general, this descriptive 
name stands for a class of situations in which EM fields 
propagate in the atmosphere from a source to a target and 
back to a receiver (detector) in the vicinity of or at the 
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same location of the source. Typical situations which belong 
to this class are optical radars, communication links (ground- 
satellite/mirror-ground) , remote monitoring of atmospheric 
pollutants, and others. 

The class of point-ahead problems such as pointer-tracker 
(with fast moving targets) , power transmission from space to 
earth (which follow a beacon signal from earth to space) 
shares a major question with the previous class, that is to 
say, the correlation between the two parts of the path. 

In this paper we are concerned with scintillation caused 
by propagation along a folded path when the target is a 
specular reflector (either a plane mirror or retroref lector ) . 
Transition cases in which the two parts of the path are par- 
tially correlated are included. To the best of our knowledge, 
only a smal], number of works in the open literature address 
this problem. We will compare the results of these works 
with ours at the related parts of our paper. 

In the second part of this^‘paper (discussion) we describe 
some accepted formulas, physical pictures and approximations, 
which describe the atmospheric turbulence and the propagation 
of EM fields in it. In the third part we develop an expression 
for the correlation between the two parts of the path. In 
the fourth we give the expressions for the scintillation 
strength (the normalized variance of the log-intensity) for 
a plane mirror and retroref lector targets. The fifth part 
describes the experimental part of our work (system descrip- 
tion and results), and a summary constitutes the sixth part. 
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II. BACKGROUND 



A. THE SCALAR EQUATION OF PROPAGATION 

The effect of the atmosphere on the EM fields propagation 
manifests itself through the Maxwell equations. The atmos- 
phere is considered to have no free charges (hence zero con- 
ductivity and unit magnetic permeability) . It is also 
assumed, based on experimental findings [8], that the time 
scale of changes of the atmospheric index of refraction is 
much longer than the time scale of the EM field fluctuations 
in the visible and near visible band. This assumption allows 
us to suppress the time dependent part of the EM field in 
Maxwell's equations. Based on the above conditions, we find 
that the vector equation for the EM field propagation is 

V^E(r) +k^n^(r)E(r) = V (V -E (r) ) (II-l) 



where : 

k = (jj/c is the EM field wave number 
n^ (r) = £ (r) 

£ = the dielectric constant 

n = index of refraction. 

It was found analytically [11,12] and confirmed in experi- 
ments [13] , that the depolarization of EM fields in the 
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forward direction due to atmospheric turbulence (variations 
in n(r)) is negligible. Thus the term, on the right hand 
side of Eq. (II-l) which represents the depolarization 
(through mixed derivatives) can be neglected [8] , and we get 
the scalar equation of propagation 



B. ATMOSPHERIC TURBULENCE 

As we can see from Eq. (II-2) , the effect of the atmos- 
phere on the EM field propagation is caused by the fluctua- 
tions of the index of refraction. The index of refraction is 
a very complicated function of the atmosphere's constitutents 
(aerosols, water vapour , etc.), the absorption coefficients 
[14] (continuous and discrete), the temperature and the pres- 
sure. Assuming the absorption is negligible (the wavelength 
does not match any absorption line or band) , the humidity 
gradient is negligible, and the aerosols' contribution to 
forward scattering is small, then the index of refraction 
turns out to be a relatively simple real function of the 
temperature and the pressure [15] 



V^E(r) + k^n^(r)E(r) 



0 



(II-2) 



n (r) 




where: 



P 



total pressure in millibars 
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T = temperature in degrees Kelvin 
A = wavelength in micrometers 

Even though the expression for the index of refraction is 
a simple one, it is impossible to get an exact solution for 
the scalar equation because we cannot specify the values of 
the temperature and pressure along the path. The only way 
we can solve this equation is through the use of their sta- 
tistical properties. 

This approach leads us to express the index of refraction 
in the following form 

n(r) = Uq + nj^(r) (II-4) 

where n^ is the average value of n (r) along the path and 
nj^(r) is the spatial fluctuation. n^ is found by replacing 
P (r) , T (r ) in Eq. (II-4) by Pq, Tq respectively, and it 
obviously cannot cause scintillation. In almost all works 
about the effect of the atmospheric turbulence n^ is taken 
to be 1. 

Since the pressure fluctuations (in the microscale) are 
relatively small [16], the fluctuations of (r) are governed 
by temperature fluctuations. 

Due to their small order of magnitude (less than one 
degree) , the main processes through which temperature inhomo- 
geneities in the atmosphere dissipate, are convection on a 
large scale and molecular dissipation on a very short scale. 
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Therefore, temperature inhomogeneities (fluctuations in 
space) which have scale size larger than a certain critical 
value, are closely related to the velocity field within the 
atmosphere . 

Using this relationship between these fields and the 
Kolmogorov model of the velocity field, Tatarski [1] gave a 
statistical description of the index of refraction fluctua- 
tions in the atmosphere. 

The Kolmogorov [17,18] model is based on the following 
assumptions ; 

1. The fluid is incompressible. 

2. The fluid is in a steady state. 

3. The fluid is getting a constant power per unit mass 
from an outside source. 

4. The Reynolds humber of the flow (onset by the outside 
source) is much greater than the critical Reynolds 
number of the air. 

5. Kinetic energy dissipation into heat can be done 
only on a very short length scale. 

Since the Reynolds number of the original laminar flow 
is larger than the critical value for air, once the flow is 
disturbed, it will turn into a turbulent flow with vortices 
whose length scale is limited only by the dimensions of the 
flow. This scale is called the outer scale of the turbulence 
(Lq) , and in the surface layer it is usually taken to be the 
height of the point above the ground. 
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If the Reynolds number of the flow within this vortex is 
larger than the critical Reynolds number, the vortex will 
break down into smaller vortices. Since the Reynolds number 
is proportional to the scale size of the flow, the formation 
of smaller vortices reduces the local Reynolds number within 
each vortex until it is smaller than the critical value, and 
the flow is stabilized. The smallest vortex size is called 
the "inner scale of the turbulence" (^q) • 

In the literature about EM fields in the atmosphere, 
these vortices are sometimes called eddies. 

C. ENSEMBLE AVERAGING, STRUCTURE FUNCTION AND TURBULENCE 
SPECTRUM 

The averaging process through which we get the moments 
of the EM field should be done over many points. Each solu- 
tion of the scalar equation represents a different n^(r) 
field (the ensemble average). Since the n^(r) field has a 
relatively short correlation length, averaging over the 
distribution of n^(r) along a propagation path can replace 
the ensemble average. In the literature the ensemble aver- 
age is denoted by < >. We will use the same notation. 

Assuming the n^ (r) field is isotropic and stationary, 
we can define the following functions. 

1. The autocovariance (sometimes called correlation 
or covariance) function 

®n^-l'-2^ “ ^ (II-5) 
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2 . 



The structure function 



Dn^-l'-2^ = < [n^ (r^) -n^ (£ 2 ) ] ^> (II-6) 

3 . The spectrum function 4>^ (K) 



B„(£i,r2) 


\ 

+<» iK(r -r^) 

= / / /4>n(K)e ^ dK 

— 00 



Since we assumed isotropic turbulence, the last 
equation can be written as 

00 

B^(p) = 4 tt J (j)^(K) dK (II-7) 

In the last four equations: 



^1 = 


spatial vectors 


p 


l£l - £ 2 ! 


K 


spatial frequency vector (of the n, (r) 
field) ^ 


K 


|K| 



Using the Kolmogorov model, the spectrum function is 
given by 

(J)j^(K) = .033 C^k"^^'^^ for 2 tt/Lq < K < 2-n/l^ (II-8) 
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-r S 



2 

This interval is called the inertial range. 
represents the strength of the turbulence and it is 
proportional to the variance of n^ at a point r . 

In most early works, this spectrum was used (by way of 
approximation) in the range 0 £ K £ This procedure as 
pointed out by Yura [19] , can cause serious deviation in the 
outcome of some calculations. Thus a somewhat different 
spectrum is now being used, namely the modified Von Karman 
spectrum. 



♦„(K) = 



2 -‘“q)' 

.033C e 
n 



(II-9) 



Hill and Clifford [20] showed that a high K values, 
the is not a decreasing function, but rather it in- 

creases to a local maximiim and then goes to zero as the 
modified Von Karman spectrum. 



D. APPROXIMATE FORMS OF THE SCALAR EQUATION 

The most commonly used approximations for the scalar 
equation are the Rytov (small or smooth perturbation) approxi- 
mation, and the parabolic equation approximation. 

The Rytov approximation is based on the perturbation 
theory and has one basic assumption: 



jVlnEj^l << [VlnEpI 



(II-IO) 
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This restriction is equivalent to the assumption that the beam 
can be scattered by turbulence eddies only once, which means 
that it restricts the range of applicability of this approxi- 
mation to weak turbulence levels. 

Starting with the scalar wave equation, and replacing n 
by 1+n^, and E by exp (T) , we get 

V^T + (VT)^ + k^(l + n^)^ = 0 

2 

Taking T = Tq + T^ and n^^ << 2n^, we get 

V^Tq + (VTq) ^ + k^ + V^T^ + (VT^) ^ + 2 (VTq) (VT^) + 2k^n^ = 0 



We can see that condition (II-IO) is equivalent to 



VT^l << |VTq 



and that if e “ vacuum solution of the scalar 

equation, the last equation can be written as 



V^T^ + 2VTq-VTj_ + 2k^n^ = 0 



(II-ll) 



which can be solved (VTq is known) , for the moments of T^^. 

Experimentally it was found [7] that the Rytov approxima- 
tion is valid as long as 



2 

a 

X 



.124 k 



7/6 



^ 11/6 




.3 



( 11 - 12 ) 
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where ; 



k = wave number, 

L = distance of propagation, and 

2 

C = the turbulence strength, 

n 

The Markov approximation is based on the assumption that 
(for an EM field which propagates along the x axis) 

1^1 << (11-13) 

3 X 

where 



u = u(x,y,z) 



E (x , y , 2 ) 
exp (ikx) 



(11-14) 



As Tatarski showed [2], the condition (11-13) is equivalent 
to the condition 1 << which is almost always met in the 

atmosphere (£q '• 1mm) for the visible band. 

Inserting (11-14) in the scalar equation and using condi- 
tion (11-13) gives the following equation, which is called 
the parabolic equation 

2ik |^u(x,r^) -i- 7^ u(x,£q) + 2k^n^ (x, r^) u (x, r^) = 0 (11-15) 

where r^ is a vector in the y,z plane. 

The physical meaning of condition (11-13) is that the 
scattering due to the eddies, is primarily in the forward 
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direction [8]. Eq. (11-15) can be solved by using the Markov 
approximation [2] or by path integral. Dashen, who was the 
first to show that the parabolic equation could be solved 
by the path integral technique, argues [10] that the con- 
straints imposed on the solution of the parabolic equation by 
both techniques are the same, and identical to the one imposed 
by the parabolic equation. 

E. PATH INTEGRALS AND THE SOLUTION OF THE PARABOLIC EQUATION 

In the forties, the path integral technique was developed 

R.F. Feynman in an effort to use the principle of least 
action as the basic principle of Quantum Mechanics. A good 
description of the technique and its application in Quantum 
Mechanics and a few other fields of physics can be found in a 
textbook by Feynman and Hibbs [21] . The use of the technique 
in solving the parabolic equation is described in a report 
by D.L. Palmer [22] . 

The use of the path integral technique in solving differ- 
ential equations can be demonstrated by the following example 
(free particle in a one dimension motion) . 

Applying the Shrodinger equation to this case gives us 

2 2 

-iTi (+ boundary conditions) (11-16) 

dt. 2m 

where \p is the probability amplitude function. 

We know that for a free particle in a one dimension 
motion. 
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and the action integral is given by 



S 




L dt 



(11-17) 



Starting from Feynman's basic assumption, namely that 
the probability amplitude assigned to the ith path (connecting 
points a,b) is 



b 

i);^(a,b) = / exp [is (a,b) /h] d£^ 

a 

where the integral is being carried out along the ith path. 

Since the total probability amplitude is the sum of the 
contributions from all possible paths, it is given by 

i|;(a,b) = /d (paths) exp [is (a, b) /h] (11-18) 

Comparing Eqs, (11-16) and (11-18), we can see that the 
solution of the differential equation (11-16) can be found 
by the use of the path integral technique, providing we can 
come up with the kernel of the differential equation, which 
in our example was obvious. 

In the case of the parabolic equation, the procedure is 
more elaborate, and since this is not the main topic of our 
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work, we would rather merely quote the result in Section 
(III.C) of this work. 

F. SCINTILLATION AND SATURATION 

Viewing a source of constant output through the atmos- 
phere, we see that the intensity varies with time (scintil- 
lates) , just as in the case of stars twinkling. 

Usually the scintillation strength is measured by the 
variance of the normalized intensity 

2 

Oj = <(I/Iq)^> - <I/Iq>^ = (11-19) 

where : 



Iq = <I>. 

Since it was found that in weak turbulence, the log ampli- 
tude of the EM field is a random variable with normal distri- 
bution [2] , the log of the normalized intensity is distributed 
normally in this region. Hence it is useful to work with 
the variable I = ln(I/lQ) and to define the turbulence strength 
by 



2 



<l^> 



- <£>■ 



( 11 - 20 ) 



Defining AI = I- Iq, we can see that if AI << Iq, since 
<AI> = 0, then 
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IhU/Iq) 



ln[ (AI + Io)/Io^ 





Thus Eq. (11-20) can be written as 




<(AI/Iq)^> - <AI/Iq>^ 



2 

<AI > 



( 11 - 21 ) 



There are several models by which the scintillation can 
be explained. The simplest model which was originally used 
by Tatarski [1] , describes the eddies in the atmosphere as 
spherical lenses whose index of refraction differs from their 



It was believed that this picture gives good approxima- 
tion only in the weak turbulence region (not more than one 
scattering event) , until Yura [23] showed that it can be 
applied to the strong turbulence region also, if we incor- 
porate in this picture the notion of the transverse coherence 
length (^q) and its effect on the focusing ability of the 
eddies . Due to focusing and defocusing of the beam while 
propagating through the atmosphere, the transverse coherence 
length of it is gradually reduced and therefore, the ability 
of the eddies to focus the beam further, is reduced too. 

Hence the growth of the scintillation strength is saturated 
and it should reach a maximxam. 



surroundings by An = n^. These lenses focus or defocus 
(depending on the n^^ values) the rays onto the target plane. 
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other models which are used, either in detailed calcula- 
tions or in qualitative explanations, are the phase screen 
model [6] and the scattering model [24] . The first describes 
the effect of the atmosphere on the beam as a distortion of 
the phase front due to different optical paths in the turbu- 
lence. In turn, the distorted phase front will produce an 
interference pattern in the target plane. The second views 
the EM fields at the target plane as the sura of fields which 
were scattered at small angles by off-axis eddies and the direct 
(unscattered) field. 

Since in all models, the scintillation strength depends 
on the turbulence strength (density and effectiveness of the 
eddies) and the total number of focusing or scattering events 
along the propagation path, it is found very frequently in 
the literature that the integrated turbulence strength (the 
effective turbulence strength) is given by the Rytov param- 
eter which was defined in Eq. (11-12) . Sometimes it is 
denoted by B or 

It was found experimentally [7] that in the weak turbu- 
lence region (also called the Rytov region) , in which the 

intensity distribution is log normal, goes up to 0 . 3 . As 
2 

grows beyond this value, the intensity distribution is 

close to the K distribution [25] , the scintillation strength 

2 

reaches a maximum [26] and starts to decrease. As a grows 
further, the intensity distribution is best described by the 
exponential distribution [10,27] and normalized intensity 
variance approaches 1. 
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III. STATISTICAL DEPENDENCE OF THE TWO 
PARTS OF THE FOLDED PATH 

A. PREVIOUS WORKS 

As we mentioned in the first chapter of this work, the 
question of whether the two parts of the folded path are 
statistically independent, is a basic one when we try to 
describe the statistical behaviour of the EM field which 
propagates along it. Once we can assxime statistical indepen- 
dence, we can treat the two parts separately, applying the 
well known results for the one-way propagation to each of 
them. Since our main object is the scintillation along a 
folded path, with a special emphasis on short distances be- 
tween the source and the detector, the first subject of our 
analysis is to try to define a criterion for the statistical 
dependence between the two parts. 

Among the works that we were able to find in the litera- 
ture, which analyze the folded path question, some like the 
theoretical work by Yura [28] and the experimental work by 
Ben-Simon et al., [29], assume statistical independence of 
the two parts. Yura gives the condition for independence as 



where d is the distance source-detector and D , D, are the 

s d 

aperture diameters of the source and the detector respectively. 
Since in Yura's paper the target is diffusive (volume of 
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air molecules) and the reflected beam is totally incoherent, 
we cannot apply this result to our case (specular target) . 

In other works [30,31,32] the question of the indepen- 
dence of the two parts is related to the wavefront tilt or 
angle of arrival (AOA) fluctuations which depends on the large 
scale eddies. Since we are concerned with scintillation 
which is caused primarily by the small scale eddies, we can- 
not use their results. 

Before we continue our analysis, we would like to mention 
without going into detailed explanation, that the usual cri- 
terion for independence between two diverging paths (their 
end points should be separated by a few isoplanatic patches) , 
cannot be applied in our case, due to the fact that in a folded 
path case certain parts of it (near the target) might be very 
close and certainly correlated to some degree. As we will 
see later, this correlation affects the scintillation and 
we cannot ignore it. 

B. THE EXPRESSION FOR THE PHASE VARIANCE 

We start our analysis with the effect of the atmospheric 
turbulence on the phase of the propagating EM field. 

The phase fluctuations (S) relative to the vacuum value, 
due to the non-uniform medium, are given by the integral 
along the optical path. 

L 

S = k / n, (x)dx (III-l) 

0 
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0 



0 



since <nj^(x)> = 



<S> 



L 

<k / n, (x) dx> 
0 



L 

k / <n, (x) >dx = 0 

0 



The variance of the phase fluctuation is 



2 2 2 2 
= <S > - <S> = <S > 

O 



and 



<S^> 



= < 



(k / 



n^ (x) dx) 



(k / 



nj^ (x ' ) dx ' ) > 



L L 

/ dx / dx'<n, (x)n, (x')> 
0 0 ^ 



<S 



L L 

/ dx / dx'B (x,x') 
0 0 ^ 



(III-2) 



Since in the Kolmogorov model, B^(x,x') = where 

y = |x-x'| and it is the same function regardless of the 

2 

origin of y, <S > can be approximated by 



<S 



2 



> 



2 

k L / B (y) dy 

-L/2 ^ 



(III-3) 
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provided L >> y where y is the range at which 

(y) << 3^ (0) . 

n n 

As we can see from Eq . (III-2) , the variance of the phase 
is related directly to the field of the index of refraction. 
This is the reason we chose to utilize the phase variance in 
the criterion for statistical dependence/independence between 
the two parts of the folded path. 

In this part of the work we chose to work with an infinite 
plane mirror as a folding target for two reasons: 

1. It is easy to compare the case of a folded path (with 
a plane mirror) when the angle between the two parts 
of the path is 26, to the case of a straight path 
(29 = 7 t) (see Fig. 1) . There is no discontinuity of 
the rays and there are no singular points at the target, 
as in the case of the corner cube. 

2. When an infinite, or at least very large mirror is used, 
we do not have to take into account the possibility 
that part of the beam will not be reflected. Thus at 
this stage, we can neglect the effect of beam wander. 

C. THE RELATIONSHIPS BETWEEN THE MEANS OF THE EM FIELD IN 

VACUUM, THE EM FIELD IN A NON-UNIFORM MEDIUM, AND VALUES 

OF THE PHASE VARIANCES 

Following Dashen [10] , we can write an expression for £ 
(where E = Re(e) at a point (L,r) due to a point source at 
(0 ,£q). This analytical signal is the solution of the para- 
bolic equation [22] and it is given by 



28 






'I 

i 






f4l 



'I 



(III-4) 



£(0,rQ,L,r) = /d ( paths) 



X exp{ik / [^(r ' (x) ) (r+e x) ] dx} 

0 IX 



where; 



r' (x) 



dr(x) 

dx 



e 

— X 



= unit vector along x axis. 



From Eq. (III-4) we can get the expression for the analytical 
signal of the vacuum field by setting n^ = 0. Thus we get 



£ (Of r Q f L f r ) 



i 

2k 



/d ( paths) 



X 




/ [r'(x)]^dx} 

0 



(III-5) 



Taking the ensemble average of Eq. (III-4) we get 



<£> = /d (paths) > 



(III-6) 



exp{-^ / [r ’ (x) ] “^dx - ik / n, [r (x) +e x] dx} 
^0 0 ^ 



= ^ /d (paths) X exp{^ / [r'(x)]^dx} 



<exp{ik / n^ (r(x)+e x)dx}> 

0 1 X 
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Using Eq. (III-l) we can write 



r P 

<exp{ik / n^ (r (x) +e x) dx} > = <exp(S, )> 

0 ^ - 
P P P 

where = is and S is the phase fluctuation along the path 
P. Since it is generally accepted [33,34] that the phase 
fluctuation in the atmosphere is a Gaussian random variable 

P 

with zero mean, S^ is also Gaussian with zero mean, and 

( *) 

therefore it satisfies the following formula 
<e ■^> = exp [ = ] 



Thus 



(*) 

For Gaussian random variable a, 
standard deviation a 



with zero mean and with 



<e°‘> 



/ e'^f(a,0)da = 



/27ra 



/ e°‘exp[- -^^] da 



2o 



1 °° o 2 2 

i /■ r 2ao -a 1 j 

J exp[ 2 

/27Ta -oo 2a 



/ 2 “> , 2,2 
, ex 2 ^ _ /2 l j ej;p,. (q7° L ] 

/27ra -°o 2a 



da 



= exp[a^/2] f f(a,a^)da = exp[a^/2] 

— oo 

where f(a,x) is the normal distribution function for mean = x. 
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P 2 

r <(S^)^>T 

exp [- 2 " ^ 



Using Eq. (III-2) we get 



<exp{ik / n, (r(x)+e x)dx}> 
0 ^ 



2 L L 

= exp{- ^ / / B [[(r(x)-r(x')^+(x-x')‘‘ 

^00 ^ 



Assuming that the optical paths in our case 
main contributors to the path integral, satisfy 

|r(x) - r(x')| << |x - x'l 

Eq. (III-7) can be approximated by 



L 

<exp{ik / n, (r(x)+e x)dx}> 
0 ^ 



, 2 L L 

= exp{- / / B (x,x')dxdx'} 



Inserting Eqs. (III-5) , (III-9) , (III-2) in Eq. 

get 

<S^> 

<e (0 ,r^ ,L, r) > = (0 ,r^ , L, r) exp [- — 



(III-7) 



] ^/^]dxdx’ } 



which are the 



(III-8) 



(III-9) 



(III-6) , we 



(III-IO) 
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where <S^> is the variance of the phase fluctuations, due to 
propagation along a path of length L. 

From Eq. (III-IO) we see that the effect of the turbulent 
atmosphere on the average of a propagating EM field, can be 
expressed by the variance of the phase fluctuations along the 
geometrical path. 

It is obvious that in vacuum, apart from a phase shift 
due to the reflection, the fields E_ and E_t as described in 
Fig. 1, are the same. We can ask ourselves what will be the 
relationship between and Since in the case of 

<£ 2 l^ there is no folding at all, it is obvious that in this 
case there is no statistical dependence (relating to the 
folding). Hence, the closer <£p> to less depen- 

dence (statistical correlation between the two parts) , we 
have . 

Using Eq. (III-IO) and the observation that up to a phase 
constant come to the conclusion that the statis- 

tical dependence between the two parts of the folded path, 

2 2 

can be determined by the relationship between <S^> and 

the closer their values are, the less the two parts of the 

folded path are correlated. Therefore, in the next steps we 

2 2 

will calculate and Sp (we omitted the angular brackets 

for the sake of convenience) , compare their values, and formu- 
late a criterion for statistical independence. 

Before we go to the next step of the analysis, we have 
to evaluate the effect of the constraint which was formulated 
in Eq. (III-8) , on the applicability of our approach. 
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Figure 1, The reflected (folded path) and transmitted 

EM fields at distance 2L from a point source. 
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since the point (x,r(x)) is a point on the optical path 
through the turbulent medium (which is given by the normals 
to the wavefronts) , the constraint (III-8) is equivalent to 



dr (x) 



which means that the normals to the wavefronts are almost 
parallel to the x axis, which in turn means that only small 
angles of scattering are allowed. This constraint is the 
same as the one for the utilization of the parabolic equa- 
tion [8,10] (Eq. (11-13)), and both are met for almost all 
conceivable conditions. 

D. THE VALUE OF IN TERMS OF S^ , FOR THE VON KARMAN 

SPECTRUM ^ 

From Eq. (III-2) we can write 



Working on both integrals and using the transformation 




/ dx / dx'B(lx-x'l) 
0 0 ^ 



and 




, 2L 2L 

k / dx / dx'B(lx-x'l) 
0 0 ^ 



2L 



z 



X - X 



y + z 
2 



X 
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y = X + X ' 



y - z 

2 



x' 



J( 



X,x' 

y,z 



1 

2 





0 

/ B ( I z I ) dy + / dz 

L. 



2l^+z 

/ B^(lzl)dy} 

-z 




B ( I z I ) (2L.-2Z) + 
n ' ' 1 



0 

/ dz ( I z I ) (2L^ + 2z) 



where L. is either L or 2L. 

1 

Taking z ^ -z in the second term, we get 



S 



2 

L. 





(L . - z) B (z) dz 
1 n 



Inserting the explicit forms of B^(x) [2] and the Von Karman 

spectrum for 



B (x) 
n 



= 4it 



0 



/ c()^(K,x) K^dK 



(J)^(K,x) 



.033 C^(x) 



Assuming (x) = constant and taking y = KLq, 



we get 
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(L./L|3-2) 



L. 

1 



A / 

0 






dz 



f dy y Sin(zy) 

o' 



where; 



, -3 , 2 ^2 ^ 8/3 

A = .0332 TTk C L 

n 0 



and 



z = x/Lq . 



Using formula 2.3(11) of Ref. [35] we get 



L./L 

. 1 ^ 1/3 

A / dz (L,/L^-z) (z) 



r 0 



'1/3 



where : 



.033 2^/^ LQ^^/rdl/e) 



and K, /^(z) = Modified Bessel function of the third 

^ kind. 

Solving the integral in Eq. (III-ll) (detailed solution 
is given in Appendix A) , gives for (L^/Lq) >> 1 

, A n^/^r(5/6) 

% ' ^2/3 (III-12) 

which is the same as the expression which is given by Tatarski 

2 

[2]. Thus we can see that for large (L./L„), Sx . is approxi- 

1 0 

mately linear in and 
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E. THE EXPRESSION FOR sj 

r 

The geometry of the folded path and the notations which 
we use in this section are described in Figure 2. We assume 
that the mirror is large enough so that ; 

1. Diffraction effects can be neglected. 

2. There are no fluctuations in the amount of energy being 
intercepted by it. 

We also assume that the field at the point D (E ) can be 

r 

found by the path integral technique, using SMD as the geo- 
metrical optics path. It should be pointed out that since 
the calculation of the E value by the path integral assumes 

r 

dr (x) 

I — — I << 1 (scattering at small angles) 

the application of the technique to this type of path might 
be questionable, for the reason that it includes one point (M) 
at which the above condition is certainly not satisfied. 

Our justification for the use of the path integral tech- 
nique along this type of path, is based upon two arguments; 

1. In the case of fields propagating in vacuum, we know 

that apart from a phase constant, E = E„ . Therefore, 
we can use the path integral technique (Eq. (III-5) ) 
along the path SMD by excluding the point M from the 
integration. 
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Figure 2 . The geometry and notations of the folded 

path case. M = a plane mirror; S = a point 
source; D = detector. The distances SM and 
DM are of length L. 
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2. If the results of our analysis in the two extreme 
cases ; 

a . 0 = 0 

b. 0 tt/A 

match the results which are obtained in other ways, 
the fact that at the point M, — is infinite, has 
no effect on the applicability of the path integral 
technique to our case. 

At the end of this section we will discuss the second 



As we saw in Section III.C, in order to find E we have 

r 



where the integration is along the line SMD. Thus 



argument . 



2 

to calculate S„ by 

r 




j j dxdx' B (|x-x'|) 
SMD ^ 



2 2 ^ 2 
S^ = k / dxdx 'B(|x-x'|)+k / dxdx ' B ( ] x - x ' | ) 

" 0 ^ L ^ 




L 

X J dx' B (|x-x'[) 
0 ^ 
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where integration on the interval 0,L means integration along 
SM, and on the interval L,2L means integration along MD. 

In integrals and 1^, x' is the distance of the point 
(y,x^) from the point M along the line MD, and the notation 
|x-x'| stands for the actual distance between the points, 
which is denoted by P. 

Using the transformation y = x-L; y' = x ' -L 



I 



2 



9 ^ 

k / / dydy ' ( | y - y ' | ) 



0 



I 



1 



Since ( | x - x ' | ) is symmetric in x, x', = I^, and 

therefore 




f f dxdx ' B ( I X - X ' I ) 
0 ^ 



• (III-14) 



^ L 2L 

+ k / dx / dx'B(|x-x'|) 

0 L ^ 




From Eq. (III-2) 



I 



1 




The evaluation of I 



■3 



From Fig. 2 we can see that 
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' = (L-Xj^) /Cos (29) ; dx ' = - dXj^/Cos(20) 



y = (L-Xj^)Tan(26) 



(III-16) 



x-x'l = P = [ (L-x^) ^Tari^ (29) + (x^-x)^]^^^ (III-17) 



Hence can be written as 



Cos (29 ) 



f dx / dx^ B^(P(x,Xj^)) 



2L 



Sin^ 



Applying the transformation (x,x^^) -> (P,a) and observing that 

P Sin(29+a) 



1. X = L - 



Sin (29) 



2. x^ = L - P Sin a 



Coti29) 





3x 


3x 


x,x. 


3P 


3a 


3 - ^(p,a ^ = 


3x^ 


3Xj^ 




3P 


3a 



= P Cot(29) 



4 . For 0 < a < it/2 - 9 



0 < P < 



L Sin(29) 



Sin (29+a) 



= P. 



5. For tt/2-9 < a < tt-29 



0 < P < L . P 

— — Sin a ; 
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we get 



^2 tt/2-0 



tt-28 ^2 

+ / da / PB (P)dP} 

7T/2-e 0 ^ 



Making the transformation B = •TT-(2e+a) in the second 
integral, we can see that it is equal to the first one. 



•3 = Siffw J J PB„(P)dP 



where 



L Sin(29) 

1 Sin(20+a) 

Inserting the explicit forms of the Von Karman 

(})^(K) , and assuming that <j)^(K,x) = , we get 

tt/ 2-0 ^2. 9 -n /f; 

I, = A / da / dy J dz 2 (z'^+l) '^°Sin(yz) 

0 0 0 



.033 8tt 

n 0 

Sin(20) 



where ; 



A = 



Thus 
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1 



z = LqK 



y = P/L, 



L Sin(29) 
'l LQSin(20+a) 



Using formula 2.3(11) of Ref. [35], we get 



7T/2-0 



I 3 = / da / dy ^ 1/3 



where: 



033 2^/^ 

n 0 

Sin(20) r(ll/6) 



K = Modified Bessel function of the third kind 



r = the gamma function 



Comparing the last equation with Eq. (III-12) gives 



2^/^ L, 



Sin(20)TT^'^^r (5/6)L 



tt/2-0 
/ da 

0 



J % V3<y> 



Inserting the last equation and Eq. (III-15) into Eq. 
(III-14) gives 

Sp = 2S^[1 + f(yQ,0)] 

or by using (III-13) 
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(III-18) 



S 



2 

F 



+ f(y„ 



e) ] 



where 



f(Y0/9) 




tt/2-6 

I 



0 



da 



J K3_/3(y)dy 



(III-19) 



and 



A 



2 



^0 



22/3 

(5/6) 

L Sin(29) 



yj_ = yQ/Sin( 26 +a) 



We solved Eq. (III-19) numerically and found out that for 
small angles 9/ which is the case in most of the applications, 
the f dependence on 0 is very weak, and therefore 
f (Y q/Q) - f (Y q) • 

The result of the numerical integration is described in 
Figure 3, and we can see that 



lim f(yQ,0) = 1 (when 9^0, Yq 9) 



lim f(Yr»/9) = 0 

6-^tt/4 ^ 



and therefore 
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PATHS DEPENDENCY FUNCTION 
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Figure 3. The behaviour of f(y , 6) (Eq. III-19) as a function 
of y_ (for small 6) . 
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0 




6 = tt/4 



0 = 0 



(III-20) 



2 

These are the values of S„ for the two extreme points of 

r 



the transition region (from exact folding to statistical 
independence of the two parts of the folded path) . In 
previous works on the subject of scintillation on folded 
paths, the authors assumed either exact folding or statistical 
independence of the two parts of the folded paths. Thus, as 
we pointed out in Section III-E, we can compare our results to 
those calculated by others, only at these two points. 

F. A CRITERION FOR STATISTICAL INDEPENDENCE 

We have now enough information about the behaviour of 

2 

'Sp(L,9), in order to define a practical criterion for statis- 
tical independence of the two parts of the folded path. 
Strictly speaking, the two parts are independent only for 
9 = tt/4. Yet because f{yQ,0) decreases very fast at small 
yQ and has a very long tail (approaches slowly to zero) , we 
can define the region at which we have statistical indepen- 
dence, as the region in which f(yQ,0) £ e From Figure 3 

we see that this criterion is met (for small 9) . 



yo i * 



Since 



46 



I 



L Sin(29) 

"0 = 

for small angles 9, can be written as 

where D is the distance source-detector. Thus we see that the 
statistical independence between the two parts of the folded 
path depends on the distance source-detector, expressed in 
units of the outer scale. 

As we said at the beginning of this section, one of the 
justifications for employing the path integral technique in 
our case, even though the main constraint on its utilization 
is violated at the point M, is the observation that the re- 
sults in the extreme cases (9=0; 9 = ir/4) match the results 
which are obtained by other techniques. 

For the case of 9 = tt/ 4 we got 

S^O = ,/4) = 

which is the expected result based on physical arguments and 
the results of Section III.D. 

For the case of 9 = 0 we got 

- 0) = 23^^ = 43^ = (23j2 
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When a beam propagates back and forth along the same turbu- 
lent medixom, the phase fluctuation due to each path element 
Ax is doubled. Thus comparing the random variable S in the 
case of a path of length L (S ) to the case of a folded path 

1j 

with 9=0 (S_) , we get 

r 



s 



F 



2S 



L 



From probability theory, this gives for the variances 





which fits our results. 
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IV. SCINTILLATION STRENGTH FOR THE FOLDED PATH 



A. THE RATIO BETWEEN THE SCINTILLATION STRENGTHS FOR THE 
FOIDED PATH AND THE ONE WAY PATH 

The log intensity variance for the Rytov region for 
spherical waves is given by [1] 



oliL) 



.5 c" 

n 



(IV-1) 



In the strong turbulence region (exponential distribution 
of I) , the scintillation strength is given by the variance 
of the normalized intensity [8] 



a 



2 

I 



1 + 



.69 



/ 2.2/5 

(Oj) 



(IV-2) 



Comparing Eqs. (III-12) and (IV-l) we get for the Rytov 
region 



o2(L) 



6.31 T 5/6.2 

,3/2j^5/6j.5/3 ^ ^L 



A(k,LQ)L^/^S^ 



(IV-3) 



or for a path of length 2L 



aJ(2L) = A 2^/^ L^/^ S^L 



(IV-4) 



From Eq. (III-18) we get for the folded path 
aJ(F) = A 2^^/® L^/® S^[l + f(yQ,0)] 
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We can express the ratio between the folded path scintilla- 
2 

tion strength (a^(F)) and the one way scintillation strength 
2 

(a^ (L) ) , in the weak turbulence region, by 

tn(F) 11/fi 

R = = 2-^-^/^tl + f(y.,e)] (IV-5) 

This ratio is the parameter which we are going to measure in 
the experimental part, in order to verify our model. The 
dependence of R on y^ (for small 6) is described in Figure 4. 

B. DISCUSSION AND OBSERVATIONS 

1 . The Range of R and the Saturation of Scintillation 
Strength 

As we can see from Eq. (lV-5) and the corresponding 
figure, the range of R is 

< R < ai’/s 

Therefore when we work with a system for which 9 is close to 
zero or zero, we should expect the scintillation strength to 
be about seven times greater than the scintillation strength 
of the one-way path. 

One should be aware of this large ratio especially 
while designing a system or an experiment which is based on 
the log-normal distribution. Since the region of the log- 
normal distribution, as we mentioned in Chapter II, is given 
by 
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Figure 4. The ratio (R) between the scintillation strength for 
the folded path (a?(F)) and for the one way path 
o2 (L) ) , as a function of y„. (We assume small angles 







1 




in the case of exact folding (6 = 0) , the upper limit of 
this region can be reached along relatively short distances 
and/or relatively weak turbulence levels. 

On the other hand, since saturation can be achieved easily 
(relative to the one way case) , it provides a way of carrying 
out experiments in the saturated region without the difficul- 
ties of using very long paths. 

2 . Comparison with Other Works 



of a spherical wave which propagates along a folded path 
(plane mirror) with 0 = 0 is 



Since he used the phase screen approach [6] and the Kolmogorov 
spectrum, his results are applicable only in the Rytov region 
and there is no dependence on the outer scale. 

Using the Kolmogorov spectrum and the Mellin trans- 
form (formula 6.5(15) of Ref. [35], we get from Eq. (IV- 6) 



Smith [36,37] found that the scintillation strength 




64 -rr^ k^ / dx / K(}) (K) Sin 
0 0 ^ 




a^(F,9 = 0) 
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where 



A 3 = .033(4)2^/^ r(-5/6)(-l) Cos(57t/12)tt^ 



Making the transformation 



y = x/L 



we get 



aJ(F, 9=0) = A_ / dy[y(2-y) (IV-7) 



It is well known that for the one way path of length 



L [2] 



a^(L) 



2 , 2 



= 16 TT k 



/ / 



K(J)^(K) 



c • 2 ,K' 

Sin [ — 



(L-x) 

2kL 



]dK 



Following the same steps as in the previous case, we get 






= A, 






where 



A. = .033 4 (-l) r(-5/6) Cos(5tt/ 12) k^^^ 

4 n 
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and we can see that 




Taking y = x/L we get 



aja) 




dy[y(l-y) 

0 



(IV-8) 



From (IV- 7) and (IV-8) we get 



A / [y (2-y) ] ^'^^dy 
0 

A, / [y(l-y)]^/^dy 
0 

Taking y = 2x in the upper integral and observing that the 
integrand in the lower integral is a symmetric function about 
y = 1/2 we get 



R = 



a^(F,0 = 0) 



A, 2 45/6 /"''^x(l-x) 15/64, 

0 

1/2 c/fi 

A. 2 / [x(l-x)] '°dx 

0 



2I7/6 



which is exactly the same results as we got by using the path 
integral technique. 

It is easily seen that even though the results for 
9=0 are the same (using the two techniques). Smith's paper 
cannot give the values of R for 9/0. 
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Game et al . , [38] measured the value of R (for 9=0) 

along a path of 274m at night (using a collimated beam HeNe 

laser), to be 6.668. From their paper, it is not clear 

2 2 

whether the measurements of (L) and (F) were done 
simultaneously . 

3 . The Weight Function in the Case of Exact Folding (6 = 0) 
The relative contribution to the scintillation of 
different points along the path is described by the weight 
function. 

From Eqs. (IV-7) and (IV-8) we can see that the weight 
functions W^(x) for spherical waves are: 

5/6 

For one way W (x) = [x(l-x)] ' 

Li 

For the folded path W (x) = [x(2-x)]^'^^ 

r 

The two functions are shown in Figure 5. 

From Figure 5 we can see that in the case of the folded 

path (6 = 0) , the main contribution to the scintillation comes 

from the vicinity of the mirror. In the case of the single 

path, the main contribution is from the center of the path. 

Since it is very unlikely that we can find the proper condi- 

2 

tions under which C^(x) along an atmospheric path is constant, 
we should be careful when we compare results from the two types 
of path. 



55 



1 



I 



X 




4 ^ 



0 




p 






cn 




0 






0) 




-p 






cn 


S • 


0 






(d 


x; 


0) 






o 


-p 


-p 








0 Id 


0) 






0) 


OiT3 


Q) 


• 


4 :: 


11 




X: rH 


-p 


0) 


0) 


-p 






CD x^ 


x; 




11 


G 


-p 


4J 






•H 


x: 




G 


X 




-P CJ>t5 


fd 




cn 


•H G 


G 




-p 


0) 


^ 0 


fd 


0) 


fd 


> 


i-H 




0 




(d 


x: fd 


0 


p 


cn 




-p 




G 


•H 




fd G 


ii 


0 




rH 


Q 4 0 




cn 


P 


cd 


•H 


X 




0 


0 


nd -p 




0) 


p 


•H 


0) (d 


4J 


x: 


p 


U 


^ 0 


fd 


-p 


•H 


0) 


«H 0 






£ 


x: 


0 «H 


0) 


•*> 




a. <4-1 


0 


XI 


0) 


cn 


Q) 


p 


-p x: 




fd > 


G 


fd 


-p 


4H 


•H 


0 


a 




0 T3 -M 


cn 








G fd 






0 


cn 


fd rH 


0) 


0) 




G 


0) 


x: TS 


11 


0 


^ p 


-M -H 




•H 


XI 




0 


X 


4J :s 0) 




4H 




0 


-- x: x: 




x» 


G 


X) 


4J 


P 


fd 


G ^ 


fd 


0 




4-1 


-p cn 


Oi Cn 


0) 




fd 4J 






p 




04 G 


0) 




fd 


x: 


CL) 


• 




CJi 


0) cn 


rH 


p 


•H 


IH (D 


G 




0 


0) 


cn u 


•H 


11 


4J 




a 04 cn 




0 




•H 0) 




X 


0) 


0) 


cn p 


P 




-p 


x: 




0 


4J 


0) 


E-* 


a 

X 


Cn 


<t3 T3 



in 

a; 

u 

:3 

•H 



56 



( 






; 



4 . The Effect of Path Folding on Aperture Averaging 



The discussion and the evaluation of the scintillation 
strength in previous sections assximed that we have a point 
detector. The term point detector refers to a detector with 
aperture diameter (D^) smaller than a critical length (Dq) . 

The aperture averaging can be explained qualitatively 
in the following way [7,40] . Due to atmospheric turbulence, 
the EM field at the detector plane has a certain spatial struc- 
ture. This implies that the intensity (the modulus of the EM 
field) has a spatial structure too. In order to measure the 
full effect of the scintillation utilizing direct detection 
of intensity, the detector aperture should be smaller than 
the distance along which the intensity variations are negligi- 
ble (Dq) , otherwise the detector will average different 
intensities across its aperture. (In the language of Young 

W 

[40] , the detector will filter out the high spatial frequen- 
cies of the intensity's spatial variations.) The critical 
length is given by 



Bj(Dq) = Bj(0) e 

where B^ is the covariance function of the intensity at the 
detector plane. 

It is well known [1,2] that in very weak turbulence 
Dj'Jp = (XL) 

where A is wavelength and L is distance of propagation. 



57 



I 




In the case of strong turbulence, it is not clear 



what is the exact value of [8], The difficulty in obtain- 
ing the Dq value arises from the difficulties in calculating 
the fourth moment of the EM field. Since the transverse 
correlation length of the EM field (r^) [23] is a decreasing 

function of the distance whilst I is an increasing one, at 

r 

a certain distance the field and the intensity spatial covari- > 
ance functions are dominated by r^ . Thus, beyond certain 
turbulence strength, Dq should be a function of r^ . 

As Yura [28] pointed out, the MCF of the EM field 
is proportional to the phase variance. 

It follows from the above mentioned arguments that, 
as the phase variance increases the critical length Dq 
decreases. Since we found that in the case of the folded 
path 



g2 ^ 2^^/^[l + f(yQ,0)] 

it is obvious that while working with the folded path, we 
should utilize detector apertures which are much smaller 
than those which are used in the one way case. 
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V. EXPERIMENTAL WORK 



The main object of the experimental part is to measure 

the values of the function R(yQ) and to compare them with the 

analytical function which is given in Eq. (IV-5) . 

In order to obtain the values of RCyg ) i we have to 

measure simultaneously the values of the one-way scintillation 
2 

strength and the scintillation strength for the folded 

9 

path (ol) . These measurements should be done at several source- 

r 

detector distances (assuming that is the same for all 
measurements) . 

A. THE EXPERIMENTAL SYSTEM 

The experimental system consists of two major parts: 

1) The optical system; 

2) Data reduction system. 

1 . The Optical System 

The optical system is described in Figures 6 and 7. 

The main problems which we dealt with during the 
planning and construction phases of our work were: 

a) Selecting proper sources. 

b) Making sure that the beams which propagate along the 
one-way path and the folded path see the same turbulence. 

c) Obtaining high quality optical elements with rela- 
tively large apertures (within our budget and time 
span constraints) . 
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Source 2 




Pockels Cell Modulator 



Polarizer 



Partially Transmitting 
Plane Mirror 




Figure 6 . 



The optical system for measurements with 9 7 ^ 0 . 
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Pockels Cell Modulator 




Figure 7. The optical system for 6 



0 . 
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The sources of our optical system should meet certain 
requirements : 

a) Low noise. They should have low values of intensity 
variance in order not to obscure the intensity fluc- 
tuations due to the atmospheric turbulence along 
relatively short paths. 

b) Operation in a pulsed mode or a simple way of 
modulation. (This requirement evolves from the 
way we were planning to subtract the background 
noise . ) 

c) When passing through the same turbulence, the respec- 
tive measured scintillation strengths should be highly 
correlated. 

d) Synchronized operation either in pulsed or modulated 
modes of operation. This requirement makes it possible 
to compare pairs of pulses . 

Technically it is easy to employ a GaAs laser at the 
plane mirror end (the one-way source) in a pulsed mode, and 
a mechanically chopped HeNe laser at the other end (folded 
path source) . These sources in these modes of operation, 
were used before by the Electro-Optics group of the Physics 
Department at the NPGS in previous experiments (41] . 
a. The Sources 

It was found (assximing the sources ' intensities 
are normally distributed) that the normal intensity variances 
of the sources were: 



62 



1 ) 



HeNe laser (Spectra Physics Model 155) mechanically 

2 -4 

chopped; £ 5 x 10 

2 -4 

2) 12 elements array of GaAs; £ 9 x lo 

2 -4 

3) 8x8 elements of GaAs; S^. £ 5 x lo . This laser 
matrix which is manufactured by Laser Diode Labs, was 
equipped with a Lucite parallelepiped rod mixing ele- 
ment, in order to enhance its output uniformity. This 
mixer was designed and built at the NPGS . 

The conclusions from these measurements were: 

1) In view of the expected levels of scintillation strength, 
all three lasers can be incorporated as sources in the 
optical system. 

2) The 8x8 GaAs is a better (more stable) source than the 

12 elements array GaAs. Therefore we decided to continue 
. the experiments with the 8x8 GaAs laser only. In the 
following text, the term "GaAs laser" refers to this 
laser . 

The correlation coefficient between the HeNe and 
the GaAs lasers were measured in a configuration which is 
described in Figure 8. The GaAs in a pulsed mode was triggered 
by the mechanical chopper via a pulse generator. In this way 
we were able to get interlaced series of pulses with a con- 
trolled time separation between them. For each laser we 
measured the log-intensity variance (in a way which will be 
described in Chapter V.B) and each pair of variances (one 
for each laser) was compared. For each run we calculated the 
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HeNe laser 




GaAs pulsed laser 



Figure 8. The configuration of the optical elements for 
the measurement of correlation coefficients 
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statistical correlation coefficient between the two series 
of variances' values. One run of this part of the experiment 
is described in Figure 9 . In this particular run we got the 
highest value of correlation coefficient (.973) . In the other 
runs the correlation coefficient values were much smaller. 

As a result of this measurement, we decided to 
replace the GaAs laser with a second CW HeNe laser (similar 
to the first one) , which is modulated by a Pockels cell 
modulator. (We tried to use a second mechanical chopper in- 
stead of the Pockel cell modulator but failed, because we were 
unable to synchronize it with the first one.) We found that 
the values of the two series of scintillation strengths (log- 
intensity variances) of the two HeNe lasers were highly 
correlated. A typical run with p = .996 is described in 
Figure 10 . 

The measurements of the correlation coefficients 
were carried out in the basement corridor of Spanagel Hall 
at the NPGS. (Its structure and turbulent features are des- 
cribed in Ref. 41.) The path length was 145 meters and the 
time separation between the pulses of the two lasers was 
250-350 ys. 

Based on the results of these experiments, we 
decided to continue our experiments with the two modulated 
CW HeNe lasers. 

b. Frozen Atmosphere 

As we mentioned in Chapter II of this paper, the 
atmosphere can be considered frozen for time intervals of 
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Figure 9. Correlation between the log normal intensity variance 
of the HeNe and the GaAs lasers which propagate along 
the same path. The value of the correlation 
coefficient p is .973. 
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Figure 10. Correlation between the log normal intensity variances 
of two modulated CVJ HeNe lasers which propagate 
along the same path (p = .996) 




i 




up to 1 ms. Since the round trip time of flight along the 
folded path (when the target is up to 10 km from the source) 
is less than 70 ys, the atmosphere can be taken to be frozen. 

In our experiments we employ two interlaced 
series of pulses with a time separation (At) within each 
pair of pulses. The time separation in our system is dic- 
tated by the equipment which was available to us, and is not 
less than 250 ys. 

We tested the assumption of frozen atmosphere 
by measuring the correlation coefficients of two series of 
log-intensity variances. The two series of pulses originated 
from one source by using two series of slits on the m.echani- 
cal chopper for the HeNe laser, or by the use of the double 
output signal of the pulse generator, for the GaAs laser. 

For the "double pulse" HeNe laser experiment 
with t = 270 ys, the correlation coefficient was p= .998. 
This experiment which was carried out in "the corridor" shows 
that for time separation of up to 270 ys, the atmosphere is 
frozen . 

This experiment served also as final checkup for 
the equivalency of the two channels of the data reduction 
system (see Chapter V.B) . 

c. Optical Elements 

We originally planned on doing the experiments 
with two targets; a plane mirror and a corner cube. The 
original configuration of the optical system was dictated 
by four major requirements: 
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1 . 



High quality elements. To maintain the phase relation 
between different parts of the beam. 

2. Large aperture elements. To intercept a large part 

of the beam, in order to avoid the effect of different 
transmitted/reflected intensities on the scintillation 
strength . 

3. The laser at the target end (the one-way transmitter) 
cannot be placed behind the target, because its beam 
cannot be transmitted through the corner cube. 

4. The beam-splitters in use should not cause an inter- 
ference pattern. 

The original configuration of the optical system 
is described in Figure 11. 

The beam-splitters and the plane mirror target 
are made out of large pieces of flat fused Quartz, which were 
coated with Al, to achieve the required reflectivities, at 
the NPGS facilities. They were checked for flatness by 
interferometric techniques. A full description of the testing 
and coating processes is given in Reference 42. 

The interference patterns due to the double reflec- 
tion from the two surfaces of the beam-splitters, can be 
eliminated by applying an anti-reflection coating to one of 
the surfaces. Since we were not able to use this technique 
for various reasons, we used a different approach. 

The beam-splitter at the detector end is located 
very close to the detector (5-10 cm) . Hence, by using thick 
plates (taking apart the two secondary sources) we can eliminate 
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Figure 11. The original configuration of the optical 
system for measurements with 9=0. a 
is the Brewster angle. 
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this problem when the beam has a small angle of divergence. 

This beam-splitter thickness is 1/4", it has an elliptical 
shape (4 x 6"), and was coated to give a ratio of 1:2. The 
beam-splitter at the target end is a circular (6" diameter) 
plate of 1/2" thickness. As we can see from Figure 11, its 
interference pattern (at the detector) was eliminated by 
applying the Al coating to the back surface and working with 
a polarizer and angle of incidence equal to the Brewster 
angle, to eliminate the reflection from the front surface. 

In the course of the experiments, it was found 
by observing the beams at the detector end, at night, that the 
beam which originates at is very homogeneous but the 
homogeneity of the beam that originates at suffers a great 
degradation due to its double pass through the Quartz beam- 
splitter at the target end. Since we use a back surface 
reflection for the first beam (it passes this beam-splitter 
twice too) , the only explanation for the different effects on 
the two beams, is that the inhomogeneities in the beam- 
splitter material are of larger length scale than the diameter 
of the beam near the laser output aperture. 

Because of the beam quality degradation, we changed 
the optical configuration from the one which is described in 
Figure 11, to the one which is described in Figure 6, uti- 
lizing the circular beam-splitter which has 70% reflectivity 
as the target mirror (with the coated surface facing the 
detector end) . 
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As a result of the above mentioned change, we 
abandoned our original plan to measure the R values with a 
retroref lector as a target. 

d. Other Elements 

1) The mechanical chopper is constructed out of circular 
metal plate which is rotated at about 3000 rpm by an 
electrical motor. There are 20 slits on the circum- 
ference of the metal plate with a ratio of 1:4 between 
the slit's width and the opaque regions between them. 
The chopper modulates the CW HeNe laser beam at a fre- 
quency of about 1 kHz with "pulse" width (laser on) of 
about 200 ys. An LED (GaAs) and a detector which are 
placed on opposite sides of the chopper plate, provide 
a trigger signal at the same frequency as the laser 
modulation for the whole system. Due to this arrange- 
ment, we do not have to use a stabilized chopper. 

2) The Pockels cell modulator is made by Coherent Asso- 
ciates (Model 3003) . 

3) The diameter of the aperture in front of the detector 
is 3 mm. 

2 . Data Reduction System 

The data reduction system is described in Figure 12. 
It is a modification of a system which was used by the EO 
group for measuring the scintillation of one source. This 
system and its mode of operation are described in detail in 
Refs. [41, 42] . 
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Detector 

(Silicon APD — GE Model 50HS) 



chopper 



Amplifier 
(PAR model 113) 



Demodulator 



Demodulator 



C 
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K 



Log-Convertei’ 
HP 7562 a 



Electronic 
Switch 

AD 7510D1 



ATD 



Nicolet 
Instruments 
Nic 80 



HP 9825 



Log-Converter 
HP 7562a 



■> Printer HP9871 



Figure 12. Data reduction system 
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The elements which were added by us are: 

1. The second channel (demodulator and log converter). 

2. The electronic switch. 

3. Clock. 

4. Software (for the Nic-80) . 

The electronic switch (AD 7510D1) alternately sends 
the signals from the parallel channels to the computer. 

The clock is constructed of a set of pulse generators 
which are triggered by the mechanical chopper's LED-detector 
pair. The time sequence of the signals is described in 
Figure 13 . 

The change in the software enables us to store the 
measured log-intensity values of the two series in different 
storing spaces, for later calculation of the log-intensity 
variances . 

The main steps which are carried out in each channel 
of the data reduction system, are the following: 

1. The Demodulator samples the signal intensity and later 
(between pulses) it samples the background intensity. 

It siibtracts the second from the first and sends this 
difference to the log-converter. 

2 . The log-converter takes the log of the input signal and 
sends it to the computer input via the electronic switch. 

3. The computer samples the input signal, converts it from 
analog to digital and adds one to the proper counter. 
After sampling 16,384 pulses, it calculates the variance 
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Chopper ' s LED - Detector 
pair 



Main Trigger 



Mechanically Modulated HeNe 



Pockels Cell Modulated HeNe 



Sampling line of Channel 
I signal 

Sampling line of Channel 
I background 

Output of Channel I 



Combined Output 
(After switch) 



L 

0 



— ■■I.. L 1 I - Sampling time of Nic-80 

— I 1 I 1 > 

Time scale (ms) 

.5 1.0 1.5 2.0 

Figure 13. Time sequence of the system's signals 
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assuming that the distribution (after the conversion 
to log intensity was made by the log convertor) is 
normal. The calculated value is printed on the 
assigned printer. 

Note: When we measured the sources' (lasers') normal 

intensity variances, we bypassed the log con- 
vertors, and used the average value of the 
intensity distribution to normalize it. 

After we completed the modifications, we tested the 
two parallel channels for equality by comparing two series of 
pulses which came from the same laser. The high correlation 
coefficient between the variances of the two series (p = .998) 
proves that the two channels are equal. The results of this 
experiment are described in Figure 14 . 

B. EXPERIMENTAL RESULTS AND DISCUSSION 

The final experiment was carried out along a path of 62 
meters on the lawns of the NPGS . The height of the optical 
path above the ground was 1 . 1 meters . 

Because of weather conditions, we made measurements on 
February 27th and on March 3rd and 4th, 1982. Throughout the 
experiment, the weather was fair to clear. There were light 
winds (0-3 knots) from various directions, large portions of 
the path were shaded by high trees, and most of the time the 
grass was wet. The meteorological parameters which were taken 
from the NPGS meteorological station (measured at the NPGS) 
with other relevant parameters are sximmarized in Table I. 
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TABLE I 



Information about the conditions under 
which the final experiment was done 



Distance Source R Day Time T°F RH% Weather 

Detector (cm) 



0 



8.23 3-4-82 1237-1245 56 63 Clear 



20 



40 



70 



6.19 


2-27-82 


1626-1633 


58 


80 


Scattered 

clouds 


6.12 


3-3-82 


1447-1454 


56 


67 


Scattered 

clouds 


4.35 


2-27-82 


1239-1247 


61 


68 


Clear 



105 



3.22 



2-27-82 1053-1100 58 69 



Scattered 

clouds 



During all the runs, the wind was between 0-3 knots from 
various directions. 
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Figure 15 shows the experimental results, and the function 
R(Yq/ 9) which is given in Eq. (IV-5) . 

The R value of each experimental point is an average of 
16 measurements. Each measurement is the ratio of the log- 
intensity variances of the folded path (given distance source- 
detector) to the one way. 

The y^ values of the experimental points were calculated 
by 



y 



0,i 



D./L 

1 



0 



where is the distance source-detector (for the folded path) , 
and Lq was taken as 20 cm, to given an optimal fit with the 
analytical expression. 

Even though we have only five experimental points, it 
seems that the experimental results match the analytical 
prediction quite well and confirm it. 

The outer scale value (Lq ~ 20 cm) which gives an optimal 
fit between the experimental results and the analytical pre- 
diction, seems to be relatively low, for a path which is 1.1 
meters above the groxind in an unstable surface layer. For 
a stable surface layer, this is a reasonable value [44] . 
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Figure 15. A comparison of the measured R(yQ) values (o) 
with the analytical prediction (solid line) 
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VI . SUMMARY 



Optical scintillation along folded paths with different 
angles between the two parts of the path, is a part of the 
noise in many systems which measure the intensity of a re- 
flected signal (e.g., a system which measures pollutant 
concentrations in the atmosphere via their absorption) . 

In previous works the scintillation strength along these 
types of paths was calculated by assuming either that the 
two parts of the path are statistically independent or by 
working with an "exact folding" (0 = 0) . 

In this work we utilized the path integral technique and 
the Von Karman spectrum, to formulate a criterion for the 
statistical independence of the two parts (for spherical 
waves and a plane mirror) . This criterion was given in terms 
of the distance source-detector expressed in the outer scale 
unit , namely 



YO > D/Lg > 4 

For the same case, we developed an analytical expression for 
the ratio (R) of the log-intensity variances (folded path 
to one way path) as a function of y^, in the Rytov region 

R(Yj) = + ffYij)] 



where 
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2L (L = the distance source-target) 






f(yo> = < 






0 

'o 



Yo 1 4 



= 0 (exact folding) 



and showed that this prediction gives for the end points 
(Yq = 0,Yq = 2L) the same results which were predicted in 
other papers which dealt only with these points. Our ap- 
proach gives an analytical prediction for intermediate values 
of Yq, which others do not. 

In order to verify our predictions, we designed and put 
together an experimental system by modifying an existing data- 
reduction system, and constructing the optical system. This 
system enabled us to measure simultaneously the scintillation 
strength along a one-way and a folded path, which to the 
best of our knowledge, was never done before. 

We measured R(Yq) and found that 

1) the maximum and the minimiam values of R(Yq) fit the 
analytical prediction well. 

2) the slope of the experimental value fits well the 
predicted one, for an outer scale value = 20. 

Based on these results, we think that our analytical 
expression predicts well the scintillation strength on folded 
paths . 
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APPENDIX A 



THE SOLUTION OF EQUATION (III-ll) 
The right-hand side of Eq . (III-ll) 






L. 

1 



= A, 



dz (L^/Lq-z) z 



1/3 



V 2 <"> 



is a difference of two integrals. The first one is 



^ i /^0 



.1/3, 



ll = ^l^i/^0 J 



and the second one is 






V3, 



= *1 / 



The first integral 

Using formula 7.14(6) 



II = A^(L^/Lq) 



+ H^^3(z)K_2/3(z) ] 



of Ref. [43], we get 
r (5/6) z [K^^3 (z)H_2^3 (z) 

0 



where H^(z) is the modified Struve function. 



(A-1) 
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Equations 7.2(12), 7.2(13) and 7.5(57) of Ref. [43] give 
the following relationships. 



I (z) = I (z/2) (m+r+1) 

^ m=0 



(A-2) 






(A- 3) 



where I (z) is the modified Bessel function of the first 
r 

kind. 



H (z) = I (z/2) [r (m+3/2+r)r (m+3/2) ] "^ (A-4) 

^ m=0 

For z ->■ 0 , formulae (A-2) and (A-4) can be approximated by 

I^(z) = (z/2)Vr(r+l) (A-5) 

H^(z) = (z/2)^'^^[r (3/2)r (3/2+r) ]"^ (A-6) 

Inserting (A-3) , (A-5) and (A-6) in Eq. (A-1) , and calcu- 

lating the lower limit gives 
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I 






1 






Sin ( 27 t/ 3 ) 



= tt/ 3^/^ lim{ [(z/2)'^/Vr (2/3)-(z/2) (4/3) ] 

z->0 

X (z/2)^/^[r (3/2)r (5/6) ]'^ 

+ [(z/2)~^/^/r(l/3)-(z/2)^/Vr(5/3) ] 

X [ (z/2) (3/2) r (11/6) ]"^] }z 

= tt/ 3^^^ lim{az-b z^^^ + cz^-^^ - dz^} 

z-*-0 

where a, b, c and d are constants. 

Thus for z = 0 

z[K^^3(z)H_2/3(z) + H^^3(z)K2^3(z) ] = 0 (A-7) 

For large z, formulae 7.13(7), 7.5(52), 7.5(63) and 7.13(4) 

-2 

of Ref. [43] give, up to 0(z ) , 



K (z) 
r 



(7r/2z)^/^e“^ [1 + 



r (V 2^ ) . (2z) 
r(r- |) 



^ + 0(z'^) ] 



(A-8) 
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f 










H^(z) 



/2tt 



(A-9) 



, N-1 

- I <-i) 



m=0 



"■ "‘'"V* + 0(z=^-2N-i)] 

r (r + J- m) 



Using Eqs. (A-8) and (A-9) in the upper limit of Eq. (A-1) , 
we get (z = L^/L^ >> 1) 



z [Ki /3 (z)H_ 2/3 (z) + (z) K 2^3 (z) ] 



= (TTz/2)^/^e~^{ [1 + r - j - Ji/l! (2z)~^ + 0(z"^)] 



X - r 41 Z^( 2 z)-i 0 (z- 2 )) 



/2ir 



r (-1/6) 






+ [1 + - f - jx^ y---(2z)~^ + 0(z~2)] 



X [ Q-g — [1 - r(ll^(2z)"^ + 0 (z"^)] 

r (— ) 

" IT (5/6) + 0(z ) ] ] 



Iro , rr(ll/6) r(5/6) , r(13/6) r(ll/6) 

2^^ ^r(-i/6) " r(-i/6) r(5/6) ~ r(5/6) 



+ 0(z ^)} - e ^ [az^'^^ + (smaller powers of z) ] 



where a is a constant. 
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1 



Thus for z > 10, the last term is negligible and 



z [Ki/ 3 (z) H_2/3 (z) + Hj_^3 (z)K 2^3 (z) ] 



(A-10) 



= 1 + i,rav6i . ii -1 -2 

■^ 4 'r (5/6) + u (z ) 



Inserting Eqs. (A-7) and (A-10) in Eq. (A-1) , we get for 
Li/Lq 1 10 



A^(L^/LQ)7T^/^r (5/6) 
;273 



1 ni3A)_ _ 17 -1 

4^r(5/6) 18^ 



+ 0((L^/Lq)"^)} 



(A-11) 



The second integral 

Using formula 7.14(3) of Ref. [43], we get 



^2 = A^[-z'‘''\/3(z)1 



L./L^ 
1 0 



For the upper limit we make use of Eq. (A- 8) and we get, up 
to z~^ = 2 : >> 1) 






(j) [1 



r(17/6) -1 

r(5/6) 



0(z"^) ] 



-3 

which for L./L^ = z > 10 is almost zero (less than 10 ) . 

1 0 — 

For the lower limit, by use of Eqs. (A-2) and (A-5) 
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lim (z) 

z^O ^ 



2 Sln(4V3) liml2"''Vr(-l/3)-az3/3] 



7T 2 



1/3 



Sin(47T/3) r (-1/3) 



Thus for L^/Lq ^ 10 



7T 2 



1/3 



Sin(47T/3) r (-1/3) 



(A-12) 



We can now write the solution of Eq. (III-ll) 



A^{ 



(L^/LQ)TT^/^r (5/6) 



,, l.r(13/6) 17, , , -1 

^ 4 ^!W6) T8^ ^^i/^0^ 



L. 

1 



,1/3 



+ 0(L^/Lq) ) - sin(47T/3) r (-1/3) 

. 1/2 



} 



1 ^ /. xr. , rlrr(13/6) 17, 

-2/3 (L^/Lq){ 1 + [4tp(5/5) - 18^ (A-13) 



2 IT 



1/2 



Sin(4iT/3) r (-1/3) r (5/6) +0((L^/Lq) )} 
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APPENDIX B 



NUMERICAL EVALUATION OF f(yQ,9) 



The equation for f(y^,e) is 



A„ tt/2-9 



^ ^ ^ da / y^^^K^ /^(y)dy 

^ ^0 0 0 



' 1/3 



As we showed in Appendix A, 



/ y^^^K^/ 3 (y)dy = (y) 



lim y"^^^K./3(y) 
y^O ^ 



TT 2 



1/3 



Sin(4Tt/3) r (-1/3) 



Thus 



f(yo/ 9 ) = 



^2 //2-e 

- J da[B-y/ K^/ 3 (yi)] 



where 



B = 



TT 2 



1/3 



Sin(4TT/3) r (-1/3) 



or 



f (y 



0 



6 ) 



^[B(TT/2-e) 

^0 




e) ] 



( 1 ) 



( 2 ) 



(B-1) 



(B-2) 



89 



1 



I 

1 






k 



where 



I5 (Yo'6 ) 



/ 



7T/2-0 



da 






and 



Yj_ = yQ/Sin( 28 +a) 



Taking 



6 = 26 +a 



we get 






6 ) 



20 



4/3 

/ d6y^^ K 



4/3 



(Yi ) + 



7T/2 



tt/ 2+9 

/ d6y 






From 



Yi = Yg/Sin 6 



we get 



Yn Cos 3 

dy, 5 d6 

Sin^6 



4/3‘yi> 



90 



Therefore, for the first integral in 



d6 = - 



Sin^B 



Cos B 



dYi 



and for the second integral 



Since 



where 



dB = 



Sin^B 



Cos B 



dy. 



Sin^B 
Cos B 



Sin^B 



(1-Sin^B) 



V. 



0 



/ 2 2.1/2' 

yi(yi-yQ) 



i5(yo-3) = y< 



0 



/a y \/3<y) 

^ TT — /t 72 y 

(y -y(,) 



0 



/>= y^^\/3<y* , 
^ ,2 2 , 1/2 
yo (y -yo> 



Yq = L Sin(29)/LQ (B-3-a) 

= yQ/Sin(29) = L/Lq (B-3-b) 

Yb = YQ/Sin (tt/2-9) = Yq/Cos 9 = ^ (B-3-c) 

Integrating by parts (making use of formula 7.14(3) of Ref. 
[43] gives 
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f 



f ^ y K^/3<y) 

I —T 27172 

y y(y -yo> 



^ISec"^ (y/y^) y‘'''^K^/3 (y) 



+ / Seo"'-(y/yj|)y'*^^K3^3 



Since 



Sec"^(l) = 0, 

Sec“^(y^/yQ) = 17/2-29, 

Sec"^ (yj^/y^) = 0, 

and 

y‘‘'^\/3(y) 

has a finite value. 



l5(y0/9) = -t(Tr/2-20)y^/^K4/3(y^) + 9y^/%/3(y^^) 



/ g(y)dy - / g(y)dy 

^0 ^0 



where 



g(y) = Sec"^(y/yQ)y'^'^^K^^ 3 (y) . 



i 

0 

(y) ] 



(B-4) 
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Hence 



f (Y q'S) 



^[B(tt/ 2-9) + (TT/2-0)yy\, .-(y^) (B-5) 

Yq a 4 / j a 

Ya 

+ 0Y^^^K^/3(Yb) - / g(Y)dY - / g(Y)dY 

^0 ^0 



This the exact form of f(yQ,0) . 

From the definition of y (Eq. (B-3-b) ) we can see that 

a 

in all practical situations, when the path of propagation is 

parallel to the ground but not too high, we can take y >> 1. 

a 

Using formula 7.13(7) of Ref. [43], we get for y >> 1 






(tt/2) 



1/2 5/6 

Y e - [ I 
m=0 



r(ll/6+m) -m 

mir (11/6-m) 



+ 



0( 




) ] 



= (tt/2) ^/^y^/^e"^[l + (55/72) Y~^+ 0(y"^)J 



error 
Thus 

y'^^%/ 3 (y) = (V2)^/^y^/^e"^[l + (55/72)y"^] 



By taking only the first two term.s in the brackets, the 

4/3 

in the value of y ^ K^^ 2 (y) is less than 1% for y ^ 10 . 
for y > 10 we can take 



,1/2 5/6 -y 
(2tt) ' y e 



Comparing this expression with B, we find 
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B/1000 > 



for y > 10 




coefficient (tt/2-9) , we can approximate this equation by 



N;imerical Evaluation 

The numerical evaluation of Eq. (B-6) was done in the 
Computer Center of the NPCS, on an IBM 3033. The program is 
written in FORTRAN IV and is given at the end of this 
appendix. 

The important information about the program is listed 



1) The function f(yQ/0) is evaluated with y^ as the 
independent variable. The step by which y^ is changed 
is determined by a constant (parameter) of the program 

2) The angle 9 is calculated from the equation 



3) The parameter which is given for each run, is the 
distance source-mirror in units of the outer scale. 




(B-6) 



/ g(y) dy - / g(y) dy] 

^0 



below. 



9 




L . 
1 
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\ 










4) The constants A 2 , B are constants of the program. Their 
values, which were calculated from Eqs. (III-19) and 
(B-1) , are 

A 2 = .7934 B = 1.1250 

5) The values of the modified Bessel function of the third 
kind K^(y), are calculated in different ways for the 
two regions; 

(a) 0.1 y £ 4 

(b) 4 < y 

In region (a) we used formula 7.12(21 of Ref. [43]. 

00 

K (y) = J Cosh (rt) exp [-y Cosh(t)]dt 

0 

which can be written 

K^(y) = I Cosh (rt) exp [-y Cosh(t)]dt + R(r,tQ) 



Using 



Cosh ( t) 




e-^) 



> 



1 

2 



e 



t 
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9 



J 



we get 



*30 t 

R(r,tQ) < ^ / (e^^ + e~^^)exp(- dt 

^0 

Taking 

t ^0 dx ,, j ^ , 

x = e, Xg = e, — = dt and x^ > 1, 

we get 



R(r, X q) 



< 1 /" 



dx 



X, 



for r = 4/3 ; r-1 = 1/3 and -(r+1) = -7/3 

Therefore 



00 

R(4/3,Xq) £ j / (x+1) e'^^^^dx 

^0 

1 2 2 



] , 



or 

-XoY/2 

R(4/3,X(j) £ - — (Xjj + 1 + 2/y) 

for r = 1/3 ; r-1 = -2/3 and -(r+1) = -4/3, 
Therefore 
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I 




K (h) 
r 



R(1/3,Xq) 
Since in Eq . 



< 1 

^0 

(B-6) 

^ ^0' 



-XoY/2 

e 

y 



and we calculate f(yQ,0) only for y^ ^ .1, setting 



R(r,XQ) £ 10 


gives 




^0 < 


450 for 


r = 4/3, and 




400 for 


r = 1/3. 



Thus, we conclude that' 

6 

K (y) = / Cosh (rt) exp [-y Cosh(t)]dt (B-8) 

^ 0 

+ R(r, 6) 

where 

R(r,t) £ 10 ^ = 4/3 or 1/3, and y ^ .1. 



We calculate K^(y) for .1 ^ y _< 4 by integrating Eq. 
(B-8) numerically in steps of At = 10 For y > 4, 
we use formula 7.13(7) of Ref. [43] 

= (n/2y)^''^e-y[V (y'^) 1 (B-9) 
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We take N = 10, and therefore for y > 4, the order of 
the absolute error is less than 



(tt/8) 



1/23-44-10 



10 



-8 



6) For or y^^ which are greater than 30, we use the 
approximation 



^a ' 30 

/ g(y)dy = / g(y)dy (B-10) 



For y^^ > 30 we can write 



/ g(y)dy 

^0 



30 

/ g(y) dy + R(y^) , 



where 



R(Yj_) = Sec"^ (y/y^) (y) dy 



4 / 3 , 



^1 

< (tt/2) / (y) dy 

30 



= (tt/ 2) [-y"^/2K4/3(y) ] 



30 



4/3 

For y ^ 30, y ^ 4 / 3 ^/^ ^ decreasing function. 
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Thus 



Riy^) 1 (tt/ 2) (y^ - 30)3 o'^/^K4^3(30) 



Taking 30 _< ^ 10^, 



R(30 1 y-, 1 lo"^) < (tt/2) ~ lo"”^ 



7) The numerical integration of /g(y)dy is done in steps 
of Ay = .01. 

Additional information concerning the numerical 
integration can be learnt from the program itself. 

The Program 
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CALCULATF THIRD ANO FORTH TERMS 
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